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ORIGAMI: DELINEATING HALOS USING PHASE-SPACE FOLDS 
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ABSTRACT 

We present the ORIGAMI method of identifying structures, particularly halos, in cosmological N- 
body simulations. Structure formation can be thought of as the folding of an initially flat three- 
dimensional manifold in six-dimensional phase space, origami finds the outer folds that delineate 
these structures. Halo particles are identified as those that have undergone shell-crossing along 3 
orthogonal axes, providing a dynamical definition of halo regions that is independent of density. 
ORIGAMI also identifies other morphological structures: particles that have undergone shell-crossing 
along 2, 1, or orthogonal axes correspond to filaments, walls, and voids respectively. We compare 
this method to a standard Friends-of-Fricnds halo-finding algorithm and find that ORIGAMI halos are 
somewhat larger, more diffuse, and less spherical, though the global properties of ORIGAMI halos are 
in good agreement with other modern halo-finding algorithms. 

Subject headings: dark matter - galaxies: halos - large-scale structure of Universe - methods: numer- 
ical 



1. INTRODUCTION 

Cosmological A-body simulations allow one to calcu- 
late the time-evolution of an initial density field, dis- 
cretized as a set of dark matter point particles with a 
given mass, from the distant past when the field was 
quite smooth to the present-day hierarchy of halos, fil- 
aments, walls, and voids. Identifying these structures 
remains one of the key challenges to the process of com- 
paring these simulations to observations of galaxies and 
clusters. Though some efforts have been made t o identify 
complex structures such as filaments and walls (Aragon- 
' Calvo et al.|2007[|Hahn et al.|2007l|Forerq- Romero et al.| 
2009| [Bond et al.||2010[ |Shandarin||2011l ), most of the 
focus is on identifying the dark matter nalos in which 
galaxies reside or the voids that comprise most of the 
cosmological volume. 

The two best-known ways of identifying halo s in sim- 
ulations are the Spherical Overdensity (SO, |Press fc' 

and Friends of 
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has grown q uite large (for an extensive listing, see |Knebe 
et al.||2011l and references therein), but many of these 



mo d ern methods rely at their core on the SO 
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et al. 12 006; M acieiewski et al.|2009 | |Behroozi et al.|2011| 
Elahi et al.| 2011 ), or in some otner way group parti- 
cles around density peaks ( |Eisenstein fc Hut|19 98 Stadcl 
2001] jAlIbert et al. | [20041 IN eyrinck et al .||2ti05r|Tweed 



et al. 2009p . In detail, these methods differ in terms 
ot how densities are calculated, whether they perform 
any post-processing or "unbinding" procedures, whether 
they identify sub-structures, and implementation details 
such as code parallelization. Recently the "Haloes Gone 



MAD" comparison project has tested many of these 
methods on an equal footing, finding that the differences 
for the basic halo properties in a co smological simulati on 
are well within the expected error ( |Knebe et aLp)lT] ). 

We would like to note, however, that the agreement 
between different halo-finding methods for only the most 
massive halos, or for masses defined within some ra- 
dius, is perhaps not unexpected, since one of the largest 
sources of variation between methods is the definition 
of the halo boundary or outer-edge. Often this is be- 
cause the halo boundary depends quite strongly on the 
value of a free parameter in the algorithm, such as a 
density cut-off or (in the case of fof) a lin king length 
that effectively serves as a proxy for density (Eisenstein| 



that effective l y serves as a proxy fo r density ( Eisenstem] 
fc Hut 1 1 19981 [Neyrinck et al.||2005l |Knebe ' et al.||iWllT 
Anderhalden fc Diemand||2011[ ). |Knebe et al.| ( |20lip did 



(2011fdT 

not explicitly test the agreement among halo finders of 
halo boundaries, for which we expect a criterion could be 
designed that would show surprisingly poor agreement, 
such as the maximum Cartesian halo size which w e use 
below. Going to full six-dimensional phase-space ( Die- 
mand et al.l[2006l |Maciejewski et aLl[2u09l [Knebe ei al 



201 1| Behroozi et al. 201 ip allows impressive identifica- 
tion of distinct halo subhalo cores (as well as streams), 
although even here the boundaries of h alos and subhalos 
can be ambiguous. As pointed out by |Shandarin et al.| 
(2012 1, knowledge of the initial and final conditions of po- 



sition coordinates is equivalent to knowledge of the full 
six-dimensional phase space in the final conditions for a 
classical Hamiltonian system, a fact which we exploit. 

In this paper, we present the ORIGAm]^] structure- 
finding algorithm which finds halos by testing whether 
particles have undergone shell-crossing. We set halo 
boundaries at their oute r caustic, i.e. at the out e rmost 
phase-space fold, which Zukin & Bertschinger (20101 
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have found to correspond well in an analytical model to 
a conventional density-based concept of a virial radius. 
The formation of structures in the universe has long been 
linked to the formation of ca ustics as matter piles up 



and f orms pancakes or sheets (Zel'dovich 1970 Peebles 
1980 p. 95). These caustics mark out trie boundaries ot 
multi-stream regions, i.e., locations in physical space for 
which the velocity field is multi-valued. Particles that 
have entered multi-stream regions are said to have un- 
dergone shell-crossing, and their dynamics become quite 
complicated as they settle into a bound structure flKof-| 
man et al.|1990[ [T9921 fVogelsberger e t al. 200 8] | White fc 
Vogelsberger||20091 |Shandarin| |20111 |Valageas||2011| |Vo- 
gelsberger fe White||201ip . 

To provide some intuition as to how the ORIGAMI 
method works, consider that though usually particles are 
thought of as simple blobs of mass, they can also be 
thought of as vertices of an initially regular grid (which 
is often the case for the initial conditions of TV-body sim- 
ulations; however, ORIGAMI currently does not work for 
"glass-like" initial conditions). Gravity distorts this grid, 
causing some of its cells to collapse and invert when shell- 
crossing structures form. In three-dimensional position 
space, in such shell-crossing regions, multiple cells over- 
lap at the same position. 

However, in six-dimensional phase space, these cells 
never cross, assuming a numerical-error-free simulation 
of collisionless dark matter. Instead, a three-dimensional 
manifold, or sheet, stretches and folds in six dimensions, 
forming familiar large-scale structures when the veloc- 
ity coordinates are projected out. Extrapolating to time 
zero, the grid is exactly regular in position coordinates 
and all velocities are zero, so initially this sheet is flat in 
phase space. Then in subsequent gravitational evolution, 
the sheet folds without intersecting itself in phase space. 
If there were such an intersection, then two dark-matter 
particles with different initial coordinates would have the 
same phase-space coordinates at a later time , a contra- 
diction for a H amiltonian dynamical system (Landau & 
Lifshitz 1969[). This pictu re has also recently b een ex- 
plored by |Shandarin et al.|(|2012[ ) and |Abel et aT| ( |2011| ). 
They use tetrahedral tessellatio ns on the in itial grid to 
identify shell crossings. Abel et al. (2011) use this to 



measure densities within the phase-space sheet, allowing 
for example a particularly clean visualization of the cos- 
mic web. Our framework, however, keeps track of the 
number of axes along which particles have crossed, en- 
abling structures to be classified as voids, walls, filaments 
and halos. 

Figure [l] illustrates the wrapping and stretching that 
gravity imparts to an initially flat sheet of particles in 
phase space. We plot two projections of a Lagrangian 
sheet of 256 2 particles from a full 3D 256 3 -particle cos- 
mological simulation run to redshift zero. The bottom 
sheet projects out the velocity coordinates and shows the 
familiar (x, y, z) coordinates. The top sheet offers a peek 
into velocity space, replacing the z (vertical) coordinate 
with v x /Hq, the x (horizontal) component of the veloc- 
ity scaled with the Hubble constant. Halos are visible as 
knots in the bottom sheet; in the top sheet, they show up 
as furious spikes, which if zo omed into would ideally ex- 



Implicit in the name of our algorithm is an analogy 
to origami. In cosmological structure formation, as in 
origami, the "sheet" starts out flat initially, and never 
crosses itself when viewed in phase space. But of course, 
the analogy is not complete. Cosmological sheets stretch, 
unlike origami sheets. Also, the dimensionality is differ- 
ent: the folding of an actual two-dimensional origami 
sheet occurs in three dimensions, whereas even in a two- 
dimensional universe, folding of the cosmological sheet 
occurs in four dimensions. If we flatten away the velocity 
coordinates, the situation is analogous to "flat origami," 
a restriction of origami in which the end result of the 
folding is constrained to lie flat in a plane. In the cosmo- 
logical case, dark-matter caustics correspond to creases 
in the folded structure. The field of flat origami has been 
studied mat hematically (e.g. |Hull|1994| [Lang|1996| [Hull 



2002 2006), and in the discussion section we speculate 



on some applications of those results to large-scale struc- 
ture. 

We describe the ORIGAMI method in Section [2] includ- 
ing how particles are tagged according to their morphol- 
ogy, how these particles are grouped into halos, and how 
halo properties are calculated. O ur m orphology classifi- 
cation results are given in Section [XT] In Section |3~2| we 
compare the origa mi ha lo catalog to a standard fof cat- 
alog, and in Section 3.3 we study the effects of the mass 
resolution of the simulation. Finally, we discuss partic- 
ular features and speculate on potential applications of 
origami in Section |4j and we give concluding remarks 
in Section [5l 

2. METHOD 

Most fundamentally, the ORIGAMI algorithm classifies 
particles according to their morphological structure into 
void, wall, filament, and halo particles, which we describe 
in Section [2~Tj Usually, a catalog of structures is desired, 
which requires grouping these particles. There are many 
ways to do this, and we explain our procedure, based on 



a Voronoi tessellation in Eulerian space, in Section 2.2 



hibit a spiral struc ture (e.g., Fillmore fc Goldreich|1984 



Bertschinger||1985||Widrow fe Kaiser||1993 ) 



Finally, we desc ribe our method of measuring halo prop- 
erties in Section [2.3( which will be used to compare halo 
catalogs in the following sections. 

2.1. Morphology Classification 

ORIGAMI classifies the morphologies of particles based 
on the number of orthogonal directions along which the 
Lagrangian "origami" phase-space sheet is folded. Al- 
though the origami-folding occurs in full phase space, 
for the present paper we make use only of position space, 
taking advantage of the fact that the relative positions of 
particles within folded regions have been reversed with 
respect to the initial grid. It may help for some purposes 
(such as demarcating subhalos) to include velocity infor- 
mation as well; however, as mentioned above, knowledge 
of the initial and final states of a classical Hamiltonian 
system is equivalent to knowledge of the final positions 
and velocities. The initial state for a cosmological sim- 
ulation run from an initial particle lattice is particularly 
simple and thus a natural choice of data to exploit. 

A particle's ORIGAMI morphology is determined by the 
number of axes along which particle-crossing has oc- 
curred. In one dimension, this is trivially tested by de- 
termining whether the final (Eulerian) positions of two 
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Figure 1. Distortion and folding in phase-space at the present epoch of a 128 2 -particle sheet, roughly 100 ?t — 1 Mpc on a side. This is a 
quarter of a 256 2 -particle sheet, initially a two-dimensional fiat slice through the 256 3 -particle cubic lattice. In the bottom sheet, particles 
are plotted at their familiar position coordinates, (x, y, z); in the top sheet, the z coordinate switches to v x /Ho, the x-component of the 
velocity scaled with the Hubble constant. Rows of particles are colored according to their initial (Lagrangian) y-coordinate. 



particles are out of order with respect to their initial 
(Lagrangian) positions. This means that if we index the 
particles according to their positions on the initial, reg- 
ular Lagrangian lattice, particles i and j have crossed if 
i < j but their order along that axis is swapped, i.e. their 
Eulcrian positions Xi > Xj. In two or three dimensions, 
we use this same criterion to detect crossings along rows 
and columns of the initial Lagrangian lattice. A particle 
i has been crossed along the x axis if there exists a parti- 
cle j in its initial x-oriented row such that (xi — Xj) and 
(i — j) have opposite signs, where again the indices i and 
j increase with initial s-coordinate. Note that we are 
taking advantage of the lack of substantial vorticity on 
cosmological scales; if a region in a simulation were able 
to rotate along three axes, our algorithm would detect 
that as well. 

Void, wall, filament, and halo particles are particles 
that have been crossed along 0, 1, 2, and 3 orthogonal 
axes, respectively. This number is a particle's ORIGAMI 
morphology index M. The most natural set of axes to 
use is the intrinsic Cartesian x, y, and z axes of the initial 
grid, but particle crossing may occur along other axes as 
well. In practice, using only the Cartesian axes seemed 
to detect only about half of the particle crossings, in the 
sense that sets of haloparticles that should have been 
contiguous (see Figure H below) contained non-halo par- 
ticles. We found that using three additional orthogonal 
triplets of axes filled these holes. Each triplet consists of 
one intrinsic x, y, or z axis and two 45°-diagonal axes in 
the plane perpendicular to the intrinsic axis. The mor- 



phology index M returned is the maximum M among all 
four sets of axes. 

In principle, there is a huge number of "higher-order" 
axes along which particles might cross, which are at odd 
angles with respect to the Cartesian grid. As the "order" 
of a set of axes increases, so does its minimum initial par- 
ticle separation (if a cubic initial grid is used). This in- 
creased particle separation decreases the likelihood that 
particles along these axes will interact and collapse into 
the same bound structure. Figure [2] shows a schematic 
of the initial particle separation along different choices of 
axes (here in two dimensions for clarity), including the 
Cartesian grid, axes rotated by 45° , and a "higher-order" 
set of axes. Note that the initial particle separation is 
not the same along each axis for these rotated sets of 
axes, a property which is quite rare and leads to even 
larger initial separations. 

We find that including particle-crossings detected 
along all 6 permutations of the extra set of axes shown in 
Figure [2] increases the number of halo particles by only 
3% or 4% (depending on simulation resolution), and in 
all, only about 5% of particles increase their morphology 
index M. These added particles have a negligible effect 
on the halo mass functions and a small effect on the halo 



sizes (defined in Section 2.3 1, growing the halos slightly. 
Though in principle many more higher-order axes could 
be used to find particle-crossings, these added particles 
do not change our conclusions while adding significant 
computational time to the morphology classification, so 
using only the 4 lowest-order sets of axes to determine 
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ORIGAMI morphology provides a good balance between 
strict completeness and computational efficiency. 
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Figure 2. Schematic of the initial particle separation along differ- 
ent sets of axes. The Cartesian axes (solid) and a 45° diagonally- 
rotated set of axes (dashed) have the smallest particle separations 
(shown in red). A third set of axes (dot-dashed) has the next- 
largest particle separation (shown in blue). We use the Cartesian 
grid and the three sets of 45° rotated axes (where the rotation is 
in the x-y, y-z, and x-z planes) to test for particle crossing. 

The particle-crossing detection algorithm is simple and 
very fast. For each particle, for all six axes (a;, y, z, and 
the 45°-diagonal axes in the x-y, y-z, and x-z planes), 
we march forward (not backward, since all particles are 
tested) along the initial lattice and test for any crossing 
along that axis. Denoting the number of particles in 
an initial row or column as N\y>, we test up to Nm/A 
particles in each row, i.e. proceeding along 1/4 of the 
initial lattice and stopping early if a crossing-detection 
occurs. We do not approach 1/2 of the lattice, since 
at that separation a false detection could occur as an 
artifact of periodic boundary conditions (PBC). Thus the 
number of distance tests required, in worst case (when 
a full fourth of a row is tested for each particle) is 2 x 
QNidN/4 = 37V 4 / 3 , where the total number of particles 
in the simulation is N = N^ D . The factor of 2 is from 
distance tests from PBC-corrections, and the factor of 6 
is from the total number of axes tested. 

ORIGAMI morphology can be determined at any of the 
output redshifts of the simulation, but we only use the 
positions at that redshift to test for particle-crossing. 
This means that if particles have crossed at an earlier 
redshift, we do not use that information but instead look 
for new particle-crossings. We leave the study of redshift- 
dependent ORIGAMI morphology, for example whether it 
will aid in substructure identification, to future work; 
however, we note that retaining the ORIGAMI morphol- 
ogy of the next-to-last snapshot (i.e. not allowing M to 
decrease) increases the number of z = halo particles by 
a few percent. In Section [3] we present the morphology 



results and halo catalogs using z — only. 

2.2. Grouping Halo Particles 

Once all particles have been given a morphology 
classification, we group the halo particles using a 
Voronoi/Delaunay tessellation, which provides a natural 
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ticle. A Voronoi tessellation partitions space into cells, 
such that all points inside a particle's Voronoi cell are 
closer to that particle than to any other. The Delaunay 
tessellation is the dual of the Voronoi tessellation and di- 
vides a 3-dimensional volume into a set of tetrahedra (or 
a 2-dimensional area into triangles) that connect parti- 
cles, such that the particles in two adjacent Voronoi cells 
are connected in the Delaunay tessellation. The VTFE 
density at each particle is given by <5vtfe = V /V — 1, 
where V is the particle's Voronoi cell volume and V is 
the average of V among all particles. 

We group halo particles that are connected on the De- 
launay tessellation, but to prevent long strings of halos 
from being linked, we first require that halos contain only 
one "core" or set of connected halo particles above some 
VTFE density threshold. This becomes the only param- 
eter in the ORIGAMI algorithm, which we set to 200 times 
the mean density. In this first round, any particles meet- 
ing this criterion that are Delaunay neighbors are given 
the same halo ID. We then add the lower-density halo 
particles that are connected to the halo cores by iter- 
atively adding Delaunay neighbors until all connected 
particles are associated. Finally, we group any leftover 
halo particles that are not connected to a core but are 
connected to each other on the tessellation. 

We find that the core density threshold has a small ef- 
fect on the mass functions of the grouped halos and has 
negligible effec t on the distribution of halo sizes (defined 
in Section J^3| . The effect on the mass function is shown 
in Figure [5] A smaller threshold produces fewer small 
halos and a higher threshold leads to more small halos, 
as might be expected. The differences in the cumulative 
mass functions for p cut — 150 and p cu t — 250 are within 
the 10% level compared to the value of p cut = 200 used 
throughout the paper, and the y ar e smaller than the dif- 
ference with fof. See Section [3~2l for more discussion of 
ORIGAMI vs. fof mass functions. 

The grouping procedure thus puts all of the tagged halo 
particles into individual halos. Since the tagging proce- 
dure establishes halo boundaries by identifying halo par- 
ticles, no post-processing to remove unbound particles 
is performed. In principle, there is no lower limit to the 
number of particles in an ORIGAMI halo, though it is hard 
to imagine isolated particles being classified as a halo as 
described above. However, we require that halos contain 
at least 20 particles in order to be included in our final 
catalog. 

2.3. Measuring Halo Properties 
In this section we define the halo properties used to 



compare halo catalogs in Sections 3.2 and 3.3 These 



will be used to calculate the halo properties for both 
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Figure 3. The effect of different values for the "core" density 
threshold shown as fractional residuals from the cumulative mass 
function of the p cut = 200 value used throughout the paper. Re- 
laxing this value produces fewer small halos, and increasing it pro- 
duces more small halos, though the difference is small. 



the ORIGAMI and fof halo catalogs. One of the most 
important is the definition of the halo center, since it 
affects many other calculated properties. We define the 
center to be the average position of the halo particles, 
weighted by their VTFE density. This greatly reduces 
the dependence of the location of the center (and many 
halo properties) on the low-density outer regions of non- 
spherical halos, compared to a non-weighted average. 

We determine -R200, the radius beyond which the den- 
sity drops below 200 times the critical density (p C rit)> and 
M200, the mass within this radius, by sorting the parti- 
cles according to their radius from the halo center and 
determining the maximum radius for which the density 
within this radius is greater than 200p cr it- If no radius 
meets this criterion for a specific halo, we consider i?2oo 
and M200 undefined and set their values to zero. We look 
at mass functions of M200 , which is a comm on definition 
of hal o mass found in the literature (see, e.g., |Knebe et al. 
2011 1. 



Though M200 is commonly used to define the mass of 
halos, we also look at the total mass of all halo parti- 
cles in order to include the full extent of the halo, since 
one of the key differences of ORIGAMI compared to other 
met hods is its dynamica l definition of the halo bound- 
ary. Knebe et al. (2011) found that v max , the peak of 



the rotation curve, is a stabler definition of halo mass 
across different halo finders, but this is also due largely 
to its insensitivity to halo boundaries. We use the more 
discriminating total mass of a halo in our comparisons, 
since a major motivation for our algorithm is to create 
an objective definition of the halo boundary. 

The final property we will compare is the halo size, 
which we define as the maximum diameter along the 
Cartesian x, y, and z directions. This choice is pre- 
ferred to other common measures such as -R200 because 
it captures the outer regions of the halo. The halo size 
is calculated in the following way: first we transform the 
halo particle positions to a local coordinate system to ac- 
count for PBC; we then calculate s x = max(i) - min(a;) 
for the x, y, and z coordinates; and finally we take the 
halo size to be the maximum of (s x , s y , s z ). We choose 



this definition over, for example, the maximum radius of 
the particles in the halo because it does not depend on 
the definition of the halo center. 

3. RESULTS 

We first present the morphology classification that 
ORIGA MI gi ves to particles in an A-body simulation in 
Section |3.1| We th en focus on comparing ORIGAMI to 
fof in Section |3.2[ followed by an investigation of the 
effects of resol ution on both ORIGAMI and FOF catalogs 
in Section [331 



3.1. Morphology 

In this section we use a 256 3 -particle simulation with a 
box size of 200 h~ 1 Mpc and standard ACDM cosmology 
(h = 0.73, fl M = 0.3, fl A = 0.7, n s = 1, and a 8 = 0.9). 
The results of the orig ami morphology classification de- 
scribed in Section [2~l] are shown, for a slice through the 
simulation box, in Figure [4] The middle panel shows, in 
Lagrangian coordinates, the redshift-zero ORIGAMI mor- 
phology indices M of 256 2 particles that inhabit a flat 
plane with equal z in the initial lattice. Each pixel of 
the square image corresponds to a particle. The bottom 
panel shows the particles in Eulerian coordinates, again 
colored according to M. 

The top panel of Figure |4l plotted in Lagrangia n spa ce, 
shows the VTFE density (described in Section 2.2) at 
each particle. The color scale is logarithmic, and a taint 
red contour is added at 1 + <5vtfe =61. This density, 
perhaps surprisingly low, divides halo from non-halo par- 
ticles if one (wrongly) assumes that "haloness" depends 
only on density. More precisely, the fraction of parti- 
cles exceeding this density equals the fraction of parti- 
cles that are ORiGAMi-identified halo particles. Figure [5] 
shows the four morphological components of the bottom 
panel of Figure |3J separated out for clarity. 

In these figures, the ORIGAMI identification of parti- 
cle morphologies accords with expectation, most obvi- 
ously for void and halo particles. For wall and fila- 
ment particles, the situation is harder to assess in a two- 
dimensional image, but again the classification looks rea- 
sonable. A couple of small regions appear void-like in 
this two-dimensional projection and yet are classified as 
walls; the shell-crossings producing these walls happen to 
lie in the plane of the figure. A comparison of ORIGAMI 
morphology to other morphology measures will be the 
subject of a future study. 

There is an interesting duality between the structures 
in the middle and bottom panels of Figure [4] This dual- 
ity was noticed by investi gators of the adhesion mo del of 



structure formation (e.g., Kofman et al.|1990||1992 ). H 



los in Lagrangian space are large bubbles that are qual- 
itatively like voids in Eulerian space, although Eulerian 
voids are more polyhedral, whereas Lagrangian halos are 
generally rounder. Dually, voids in Lagrangian space are 
small, as are halos in Eulerian space. The situation with 
filaments and walls is harder to see, but in Figure|4]many 
Eulerian filaments look, in Lagrangian space, like walls 
dividing bubbles (halos). 

3.2. Comparison to Friends- Of -Friends 
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Figure 4. Redshift-zero quantities measured from a 256 2 sheet of 
particles that share the same ^-coordinate in the initial-conditions 
lattice. The top two panels are shown in Lagrangian coordinates, 
in which each particle is a pixel in a 256 2 image. Top: Voronoi- 
tessellation density estimates, at each particle (see text for expla- 
nation). A red contour is drawn at 1 + 5 = 61 (sec text for details). 
Middle: ORIGAMI morphology indices M. Bottom: Here the M- 
colored particles (0-3 are shown in black, blue, yellow, and red) are 
plotted at their (x, y) Eulerian coordinates. 
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Figure 5. The four components of the bottom panel of Figure|4] 
separated into different panels. 



As p art of the "Haloes Gone MAD" comparison 
project (Knebe et al. 2011), ORIGAMI has been shown 
to be in general agreement with most of the standard 



halo-finders in use today. To look at the details of how 
ORIGAMI works, particularly in the limit of very small 
halos where methods tend to disagree the most, we cre- 
ate ORIGAMI and FOF catalogs for a 100 h~ 1 Mpc iV-body 
simulation with similar cosmology parameters as above 
(except n s = 0.93 and er 8 = 0.81) and at two differ- 
ent mass resolutions. The high-resolution simulation has 
512 3 particles and will be referred to in what follows as 
the 512 simulation. The low-resolution simulation has 
256 3 particles with initial conditions down-sampled from 
the 512 simulation, such that the only difference between 
the initial density fields of the two simulations is the res- 
olution, and it will be referred to as the 256 simulation. 

The FOF halo-finding method is one of the most widely- 
used and well-understood halo finders. It works by con- 
necting particles that are separated by a distance smaller 
than some linking length, which is a parameter that we 
set to the typical value of 0.2 times the mean inter- 
particle separation. Though popular, FOF is not without 
issues. As noted widely in the literature, this method 
of linking particles into halos causes some fraction of 
the halos to be "bridged" such that two density peaks 
connected by a filament of particles are g rouped into 



the same halo (see, e.g., Lukic et al. 2009). Addition 



ally, there is some doubt as to whether FOF captures the 
full a mount of particles involved in the collapse of the 
halo ( |Anderhalden fc Diemand 2011). This issue of the 
defini t ion of halo edges is partly what motivated |Knebe"| 
et al.| ( |2011[ ) to prefer Umax as a stable quantity by which 
to compare halo catalogs, as it largely ignores the outer 
parts of halos. 

We first compare the cumulative mass functions. There 
is good agreement between ORIGAMI and FOF in the to- 
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tal mass (Figure pb, with ORIGAMI finding more very low 
mas s ha los (especially for the 512 simulation; see Sec- 
tion |3.3[ ) and fof finding slightly more very high mass 
halos! This situation changes when considering M200 
(Figure [7]); ORIGAMI finds fe wer low mass halos, as no- 



ticed in jRnebe et al 
halo mass is used. I 



(2011), when this definition of the 
11s is because there are many more 
ORIGAMI halos that do not have a defined M200, having 
no radius for which the d ensi ty is greater than 2Q{)p cr n 
(as described in Section 2.3). This is largely due to 
ORIGAMI halos being in general more diffuse and non- 
spherical than fof halos and therefore more sensitive 
to the location of the halo center. As seen in Table [T] 
ORIGAMI starts with more halos than fof at both res- 
olutions - and even more considering halos with fewer 
than 20 particles, though we do not consider them here 
- but ORIGAMI loses a much higher percentage of its halos 
when counting only those that have a defined M 2 qo- 



ORIGAMI 512 
ORIGAMI 256 
FOF 512 
FOF 256 
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Figure 6. Cumulative distribution functions of total mass for 
ORIGAMI and FOF halos at both mass resolutions. 
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Figure 7. Cumulative distribution functions of Af200 for ORIGAMI 
a nd FOF halos at both m ass resolutions. The theore tical curves 
of |Sheth fc Tormen| ( |i999} and |Press fc Schechter| ( |1974} are plotted 
for reference and are in good agreement. 

We now turn to the distributions of halo size , de fined 
as the maximum Cartesian diameter in Section |2.3| and 
find that ORIGAMI and fof have very different distri- 



Table 1 

Number of halos for different halo-finders and simulation 
resolutions 



Halo- finder 


Sim. 


N h (> 20 particles) 


Nh (defined M200) 


ORIGAMI 


256 


33266 


21632 


FOF 


256 


32052 


29612 


ORIGAMI 


512 


269192 


124870 


FOF 


512 


178397 


156686 



butions (Figure [8]). For the 256 simulation, fof halos 
have a distribution around 0.3 /i _1 Mpc while ORIGAMI 
halos have a slightly wider distribution with a mean of 
around 0.7 /i"" 1 Mpc. Since the total mass functions are 
similar, it appears that ORIGAMI halos are more diffuse 
and extend beyond the outer boundaries of fof halos 
(which may be mi ssing these edge particles ( Anderhalden 



fe Diemand||201l]) ), while at the same time ORIGAMI ha- 
los sometimes (though not often) contain within them 
non-halo particles tagged w ith M — 2 (filament) instead 
of M = 3 (halo; see Section |2.1[ ) . Recall that no unbind- 
ing procedure has been performed on either halo catalog, 
which would potentially remove the unbound particles 
from fof halos, similar to how the interloping M = 2 
particles are ignored in the ORIGAMI method. 



0.10 

0.08 
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0.00 
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Figure 8. Distribution of halo sizes (defined as the maximal di- 
ameter among Cartesian directions) for ORIGAMI and FOF halos at 
both mass resolutions. 

The other effect causing the difference in halo sizes is 
that the very small FOF halos often contain a mix of 
M = 2 and M — 3 particles, such that fewer than 20 
M = 3 particles are grouped and therefore the group 
does not make it into the ORIGAMI halo catalog. This can 
be seen in Fig. [9j where a small region of the simulation 
containing one ORIGAMI halo and a few small fof halos is 
plotted showing fof halo particles in red, ORIGAMI halo 
particles (in which halos have at least 20 particles) in 
purple, and ORiGAMi-classified M — 3 particles in green. 
Though some M = 2 particles are close enough to their 
M = 3 neighbors to be grouped by the fof method, they 
do not meet the ORIGAMI definition of a halo particle. 

The VTFE density distributions for fof and ORIGAMI 



particles are shown in Figure 10 for the 256 simulation 



There is a more abrupt distinction between halo and non- 
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37 38 
x (Mpc/h) 

Figure 9. Particle locations in a 5 X 5 /i _1 Mpc region of the 
256 simulation, 40 h~ 1 TsApc deep in the ^-direction. Over-plotted 
as larger dots are FOF halo particles in red, M = 3 halo-classified 
particles in blue, and lastly ORIGAMI halo particles (with at least 
20 particles per halo) in purple. The small FOF halos contain some, 
but fewer than 20, M = 3 particles, and so do not become ORIGAMI 
halos. 



halo particles for the FOF distributions than there is for 
the ORIGAMI distributions, which accords with the com- 
mon conception that the FOF linking length parameter 
correspo nds to a density (b ut for a detailed study of this 
issue see 



More et al.||2011 ). It is interesting here to note 
that the double-peaked histogram of log(l + <Vtfe) i s 
a particular feature of the mass-weighted nature of the 
density calculated from the Voronoi tessellation. As mass 
resolution increases and more small-scale features are re- 



14 



solved in a sim ulation, this peak grows (see, e.g., Fig. 
in Section 3.3), whereas an Eulerian (volume- weighteaj 
measure of the density would wash out these features 



4x10 




3x10 



2x10 



1x10 



2 4 
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Figure 10. Distribution functions of the VTFE particle density 
for FOF and ORIGAMI halo and non-halo particles from the 256 sim- 
ulation. 



In Figure [TP] the VTFE density distributions are plot- 
ted for ORIGAMI particles separately based on morphol- 
ogy tag. As expected, the morphology tag is correlated 
with density, though there is much overlap, and in par- 
ticular the M = 2 particles have a long tail out to high 
densities. This tail likely corresponds to M — 2 particles 
clustered near M = 3 particles that become (small) FOF 
halos but contain too few M — 3 particles to become 
origami halos (see Fig. [9|. 



4x10 s 
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Figure 11. Distribution functions of the VTFE density for par- 
ticles tagged as M = (void, blue), M = 1 (wall, green), M = 2 
(filament, purple), and M = 3 (halo, red), for the N = 256 simula- 
tion. The distribution function for ORIGAMI halos is different than 
that for the M = 3 particles because of the requirement that halos 
contain at least 20 particles. 



3.3. Resolution Effects 

In general, we expect that as resolution increases, the 
boundaries of large halos remain fairly constant and new, 

smaller halos with lower masses are identified. Indeed, 

this is what ORIGAMI finds, as shown in Figure [12J which 
shows the ORIGAMI morphologies of the same redshift- 
zero Lagrangian sheet in the 256 and the 512 simulations. 
Recall that the 256 simulation is the same as the 512 sim- 
ulation, but with coarsened resolution on the initial grid. 
New halo regions appear in the 512 simulation, particu- 
larly in void and wall regions of the 256 simulation. 

However, the effect of resolution does not appear to 
be quite the same for the ORIG AMI and FOF halo cat- 
alogs. As noted in Section |3.2| the difference between 
the ORIGAMI and FOF cumulative total mass functions 
is much greater for the 512 simulation than for the 256 
simulation (see Figure [6]) , with ORIGAMI finding notice- 
ably more low mass halos. This means that the effect 
of increasing the simulation resolution adds many more 
ORIGAMI halos than FOF halos. This effect is not seen in 
the M 2 oo mass function (Figure [7]) because these smaller 
halos continue to have an undefined A/200 j and so the 
effect of increased resolution for both catalogs is the ex- 
tension of the M200 mass function down to lower masses. 

The reason why there are many more low mass halos 
in the 512 simulation is because a greater percentage of 
particles (56.6%) are tagged as M = 3 halo particles, 
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Figure 13. Distribution functions of the VTFE density for par- 
ticles tagged as M = (void, blue), M = 1 (wall, green), M = 2 
(filament, purple), and M = 3 (halo, red), for the N = 512 simu- 
lation. 

3.5X10 6 
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Mpc/h 



2.5x10 
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1.0x10 



5.0x10 
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Figure 14. Distribution functions of the VTFE particle density 
for FOF and ORIGAMI halo and non-halo particles from the 512 sim- 
ulation. 



Figure 12. The ORIGAMI morphologies M, shown in Lagrangian 
coordinates, of initially two-dimensional sheets of particles in the 
256 (top) and 512 (bottom) simulations. In the top panel, each 
of the 256 2 pixels corresponds to a particle. In the bottom panel, 
the 512 2 pixels are colored according to the average ORIGAMI mor- 
phologies of the two 512 2 sheets that comprise the 256 2 sheet in 
the top panel. Half-integers occur in the bottom panel because of 
this averaging. For clarity, we show a slice with a relatively small 
fraction of halo particles. 

compared to 47.7% in the 256 simulation. This is par- 
tially related to the greater percentage of high density 
particl es i n the 512 simulation in general, as show n in 
Figure 13 compared to the 256 simulation in Figure [TT) 
This increase in density also affects the FOF halo catalog, 
however it is to a much lesser degree: 47.3% of parti- 
cles are in FOF halos in the 512 simulation, compared to 
42.7% in the 256 simulation (see Figure 14 1. 



If we look again at the distribution of halo sizes in Fig- 
ure[HJ we see that not only do both ORIGAMI and FOF size 
distributions shift toward smaller halos, but additionally, 
the ORIGAMI distribution becomes wider. We interpret 
this widening of the ORIGAMI size distribution as being 



linked to the relative increase of low mass halos in the 
512 simulation compared to FOF halos (Figure [6]). The 
very small groups which had a mix of M — 2 ana M = 3 
particles in the 256 simulation, such that fewer than 20 
M = 3 particles were grouped, have a greater percent- 
age of M = 3 particles in the 512 simulation, so more of 
these small groups are counted as full halos containing 
at least 20 particles. 

As a final note regarding resolution, we mention that in 
ORIGAMI, the resolution with which structures are found 
need not correspond to the initial Lagrangian particle 
spacing, but could be a multiple of it. This could be 
useful in building a sort of hierarchical morphology tree 
from a single high-resolution simulation. Also, if there is 
much initial small-scale power in a high-resolution sim- 
ulation, producing undesirably many small ORIGAMI ha- 
los, it could be useful to increase the Lagrangian resolu- 
tion used for ORIGAMI morphology-tagging. 

4. DISCUSSION 
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The current version of the algorithm requires a De- 
launay tessellation to group the halo particles into con- 
stituent halos, but we would like to stress that ORIGAMI 
is fundamentally a particle morphology tagger. The 
tagged halo (and filament, wall, and void) particles can 
be grouped in different ways; the tessellation method 
merely was expedient and sufficient for this paper. A 
Lagrangian-space method was also tried, which grouped 
halo particles that are connected on the initial La- 
grangian grid as opposed to the final-conditions tessel- 
lation, but it was found that too many halo particles are 
connected on the Lagrangian grid which later form dis- 
tinct structures. One possibility to improve the particle- 
grouping step would be to add bookkeeping that keeps 
track of which other particles a given particle has crossed 
paths with; though it is not necessarily true that these 
particles all end up in the same bound structure, the 
information could be useful. 

Though a comparison of the ORIGAMI morphology clas- 
sification to other methods is left to future work, we are 
encouraged by the results so far: by eye, it looks as if 
ORIGAMI does very well in identifying filament, wall, and 
void particles as well as halos. The task of grouping 
these particles into individual filaments, walls, and voids, 
however, will likely prove more challenging than for ha- 
los. It may be worthwhile to use ORIGAMI in concert 
with another morphology identification method that can 
be modified to take advantage of the ORIGAMI particle 
classification. 

origami is successfully able to calculate halo catalogs 
for cosmological dark matter TV-body simulations, how- 
ever, there may be some limits to its applicability. Due 
to its nature, origami is unable to find groups in ob- 
servations because it relies on information about the ini- 
tial state of the system. Similarly, it would be difficult, 
though possible in principle, to apply origami to simu- 
lations with baryons or with irregular (such as 'glass') 
initial conditions because the particles aren't initially 
aligned on a regular grid. However, note that the idea of 
using caustics to identify the outer regions of groups has 
already been applied to measuring the masses of galaxy 



clusters (Diaferio 1999 Diaferio et al. 2005). Finally 



ORIGAMI is currently unable to distinguish subhalo par 
tides from their parent halo, though it is possible that 
using velocity information will allow ORIGAMI to identify 
substructure. 

As mentioned in the introduction, the mathematical 



field of paper origami has dev eloped quite recently ( Lang 



1996 


Hull 


1994 


2002 


2006) 



In so-called flat origami, 
the final product is constrained to lie in a plane after 
folding, likely with some regions in which many layers of 
paper overlap. A set of creases in the initial paper is only 
foldable into a flat origami design if the creases obey a 
set of laws. 

Some of these laws are irrelevant to large-scale struc- 
ture because of the different dimensionality and the in- 
homogeneous stretching that occurs, but a law that may 
be relevant is two-colorability. In two-dimensional flat 
origami, polygons in the paper must end up facing ei- 
ther up or down, so the tessellation formed by the initial 
crease pattern must be two-colorable, one color paint- 
ing the "up" polygons, and the second color painting 
the "down." In analogy, Lagrangian space can be tessel- 



lated with three-dimensional regions bordered by caus- 
tics. In this case, we speculate that regions could be 
successfully colored according to their two possible ori- 
entations or chiralities, i.e. according to whether the 
three initially right-hand-oriented axes are left- or right- 
hand oriented. This issue will be the subject of a future 
study: two-colorability would be quite special for a three- 
dimensional tessellation, for which arbitrarily many col- 
ors (not only four , as in two dim ensions) can in principle 
be required (e.g., Wilson||2002 ) . 



5. CONCLUSION 

The ORIGAMI structure-finding algorithm identifies 
particles as belonging to a halo, filament, wall, or void 
by determining whether they have crossed paths with 
their Lagrangian neighbors along 3, 2, 1, or orthogonal 
dimensions, respectively. In the present implementation, 
halo-classified particles that are connected on a Delaunay 
tessellation are grouped into individual halos, though in 
principle there are many ways in which halo particles may 
be grouped. Long strings or structures of over-connected 
halo particles are prevented by first requiring that halos 
contain at most one halo "core," or set of connected halo 
particles, that are above some density threshold. 

We have compared origami halo catalogs to fof cat- 
alogs at two different mass resolutions and found that 
the mass functions agree very well; more comp arisons 
to other ha lo-finding algorithms are given in |Knebe| 
et al. (20111. Though the mass functions largely agree, 
ORIGAMI halos are in general much larger than fof halos, 
suggesting that ORIGAMI halos are a bit more diffuse and 
spread out. Additionally, the smallest fof halos contain 
many particles that ORIGAMI classifies as belonging to 
a filament, resulting in ORIGAMI finding fewer small ha- 
los. This result stresses the difference between a density- 
based definition of a halo and the ORIGAMI method of 
detecting the occurrence of shell-crossing. 

The effect of increasing the simulation resolution is in 
general as expected, with more structures popping up at 
smaller scales. Both the ORIGAMI and fof methods find 
a larger fraction of halo particles at a higher mass resolu- 
tion, but this effect is greater for ORIGAMI. For the mass 
function, this results in ORIGAMI finding an increased 
fraction of low mass halos compared to fof. Addition- 
ally, the ORIGAMI distribution of halo sizes gets wider 
as resolution increases, while the fof distribution stays 
roughly the same shape, and both shift to smaller sizes. 
This resolution effect is not seen in the M200 mass func- 
tion, suggesting that it largely occurs at the outer edges 
of the halos and for low-density halos, both of which have 
little effect on A^oo- 

With ORIGAMI, we have defined the boundary of a halo 
as being located at the outer phase-space caustic. This 
has produced some interesting differences in the sizes and 
general characteristics of ORIGAMI and fof halos, which 
diagnostics that depend only on the central cores of halos 
tend to miss. In other algorithms, the definition of a halo 
edge usually depends strongly on a free parameter, but 
in ORIGAMI the definition of halo particles is parameter- 
free. Instead of depending solely on density, ORIGAMI 
finds halos by looking for folds in phase-space. 
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